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Abstract 


This report details the development of a new two-equation turbulence closure model based on 
the exact turbulent kinetic energy k and the variance of vorticity, The model, which is 

applicable to three dimensional flowfields, employs one set of model constants and does not 
use damping or wall functions, or geometric factors. 
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1 Introduction 


The use of Computational Fluid Dynamics, CFD as an engineering design tool 
has been accelerating in recent years primarily due to improvements in three major 
areas. The first is grid generation, where there are now a multitude of grid generation 
software packages available which can generate very complex three-dimensional grids. 
The other two areas are somewhat inter-related. Namely, computational speed and 
numerical algorithm efficiency. A combination of these two have reduced computa- 
tional run times for two-dimensional flowfields from hours to minutes and have scaled 
three-dimensional calculations from days or weeks to hours. As the above three areas 
are continually refined the cost effectiveness of CFD as a design tool should contin- 
ually increase. However, one of the major limitations of existing CFD solvers still 
remains. Namely, the inability to accurately solve complex turbulent flowfields. And, 
since the majority of flows of practical engineering interest are turbulent, this has 
become one of the most significant hindrances in the use of CFD as an engineering 
design tool. 

CFD is being applied to a wide variety of applications: jet inlet design, weather 
prediction, High Speed Civil Transport (HSCT), National Aerospace Plane (NASP), 
missile aerodynamics, combustor design, and high angle of attack airfoils, just to 
name a few. One common feature of all of the above is that an accurate prediction 
of drag, heating, and mixing rates is paramount in the design/analysis process. In 
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turn, these parameters are highly dependent on an accurate solution of the turbu- 
lent flowfield. Inaccurate turbulent solutions can cause order of magnitude variations 
in these design parameters. Therefore, the goal of the current work is to develop a 
computational model which improves the accuracy of CFD solutions when analyz- 
ing extremely complex turbulent flowfields. Concentration will be focused on flows 
which contain a combination of free shear layers, shockwaves, large adverse pressure 
gradients, and large separated regions. 

Before a discussion of the solution can begin a general description of turbulence 
is warranted. The most notable feature of a turbulent flowfield is seen by examining 
the velocity, or pressure, at a point in space, and noting that they are not constant 
in time. Rather, they exhibit very high frequency, irregular, three-dimensional fluc- 
tuations. Note that these fluctuations are three-dimensional, even for flows with a 
mean (time averaged) velocity field which is two-dimensional. The fluctuations cause 
three primary difficulties to arise when attempting to describe a turbulent flowfield. 
First, because it is a largely random process, we must rely on statistical methods 
to describe a turbulent flowfield. Second, we know that turbulent flows are highly 
rotational and three-dimensional. This inherent three-dimensionality is seen by ex- 
amining the vorticity equation and noting that the vortex stretching term (absent for 
two-dimensional flows) is required to maintain the rotationality seen in a turbulent 
flow. This implies that there can be no adequate two-dimensional approximations to 
describe a turbulent flowfield. Finally, turbulence is not a fluid property (e.g. molec- 
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ular viscosity or thermal conductivity), rather, it is a property of the flow. Therefore, 
the initial/boundary conditions will play a major role and an a priori description of 
turbulence is not possible. Note that the preceding discussion was by no means an 
attempt to completely describe the physics behind a turbulent flow. Instead, it is an 
attempt to explain the difficulties associated with the CFD solution of a turbulent 
flowfield. For a more complete description of the physics of turbulence see Tennekes 
and Lumley 1 . 

To date, there are three basic approaches to solving a turbulent flowfield. The 
first, Direct Numerical Simulation (DNS), is a time dependent, three-dimensional 
solution of the Navier-Stokes (NS) equations. In order to accurately capture the 
turbulent features of the flowfield, all of the relevant scales (time and length) must be 
resolved. An estimate for the number of grid points required 1 to capture the spatial 
scales is typically given a s Re*, and therefore, for practical Reynolds numbers on the 
order of 10 6 , this method is obviously not a feasible approach. It should be noted, 
however, that DNS is extremely useful for a limited range of low Reynolds number 
cases and, therefore, can be used as an invaluable tool in the testing/development of 
more practical turbulence models. 

The next approach, Large Eddy Simulation (LES), involves splitting the flowfield 
such that the larger scales are computed directly and the smallest scale are modeled. 
As a result LES is also a relatively expensive calculation and turbulence modeling is 
still required for the smallest scales. 
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The final method does not attempt to resolve the turbulent scales. Rather, a 
model is derived which predicts the ’effect’ of the turbulence on the mean (time aver- 
aged) flow. The ’effect’ of the turbulence on the mean flow is determined by employing 
the Reynolds Averaging procedure, whereby, each quantity of interest (velocity, pres- 
sure) is decomposed into a mean and a fluctuating component. After substituting the 
decomposed variables into the instantaneous Navier-Stokes (NS) equations and ap- 
plying a time average, two additional unknowns are produced due to the non-linearity 
of the NS equations. These unknowns, the turbulent stresses and the turbulent heat 
fluxes, represent the effect of the turbulence on the mean (time averaged) equations. 
The problem of turbulence modeling is, therefore, to devise a method of accurately 
calculating these two unknown quantities. 

There are a variety of methods which can be used to calculate the effect of the 
turbulence on the mean flow (Algebraic, One-Equation, Two-Equation, Reynolds 
Stress). We chose a two-equation turbulence closure model so that, at a minimum, the 
following criterion are met. First, the model should be ’’complete”, where complete 
refers to the fact that we only want to have to specify the initial and boundary 
conditions. This criterion requires that the closure coefficients must be constant (or 
a function of turbulent Reynolds number, R t ) for all of the flows considered. Next, 
the turbulence model must be able to solve both free shear (wakes, jets, and mixing- 
layers) and wall bounded flows. General design problems typically will include a 
combination of the above flow regimes for a given problem and it is not feasible to 
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have to use different turbulence models for each of the different regions of the flowfield. 
Finally, the turbulence model should be able to solve both two and three-dimensional 
flows and should be computationally efficient when used in either a boundary-layer 
or Navier-Stokes solver. 

Given the above requirements, and noting that existing two-equation models are 
unable to meet these criterion, we chose to develop a new turbulence closure model. 
The decision to develop a new model, rather than modify an existing model, was 
influenced primarily by the following limitations of existing two-equation models. 
First, linear k — e and k — ui models are incapable of describing wall bounded and free 
shear flows using the same set of model constants and boundary conditions 2 . Also, 
existing two-equation models are, in general, unable to predict separated flows 2,3 . 
Finally, existing models were originally developed for high R t flows, but are being 
employed in situations where R t is low through the addition of asymptotic expansions 
or ad hoc terms. Based on the above, we chose not to modify an existing model, but 
rather, to develop a new model based on the exact turbulent kinetic energy (k = 
and the variance of vorticity, or enstrophy (£ = 

The development of the current model is shown beginning with a derivation of the 
exact k-( ' equations. Then the modeling is provided along with an explanation of how 
the closure coefficients were evaluated. Once the development has been presented the 
capabilities of the k-Q model are examined, and modifications made, by considering 
a variety of flows. First, several free shear flows are considered and the k-Q results 
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are compared with the k-e model as well as with an abundance of experimental data. 
Next, a flat plate and its wake are considered. The flat plate results are compared 
with the near wall measurements of Schubauer 4 and Laufer 5 and the wake results were 
compared to the experimental data of Pot 6 and Weygandt and Mehta 7 . Additionally, 
a variety of two-dimensional airfoils are examined and the results are compared with 
the Johnson-King model 8 , the k-ui model 2 , and with experiment. Finally, the model 
is examined for a supersonic three-dimensional Cylinder-Flare and is compared with 
the experimental data of Wideman et al 9 as well as the k-e model. 
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2.1 


Turbulence Modeling 

Reynolds Averaging 


The instantaneous Navier-Stokes equations in a compressible medium are given 


as 
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with the following constitutive relations 
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(2.1) 

( 2 . 2 ) 

(2.3) 


(2.4) 

(2.5) 


Pri dxj 

where tj ix p,, Sij, and Pr ; are the viscous stress tensor, molecular viscosity, strain rate 
tensor, and laminar Prandtl number, respectively. It should be noted that in order to 


simplify the current discussion, only incompressible flows (constant p) are considered. 


For an incompressible medium the Navier-Stokes equations (2. 1-2.3) reduce to 


dui 

dxi 


= 0 



+ puj 


dui 

dxj 


dp dtp 

dx x dxj 


( 2 . 6 ) 

(2.7) 


and the constitutive (Eq. 2. 4-2. 5) relations become 
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It should be noted that the energy equation was omitted from the previous equa- 
tion set. This is due to the fact that, for an incompressible flow, the energy equation 
is completely decoupled from the conservation of mass/momentum and, therefore, 
will be excluded from the current discussion. An in depth discussion of the energy 
equation and its relation to solving turbulent flows can be found in Refs. [10, 11], 

The solution procedure for the current work does not attempt to resolve (spa- 
tially or temporally) all of the high frequency fluctuations which occur in a turbulent 
flowfield. Rather, the goal is to calculate the mean (time averaged) effect of these 
fluctuations on the mean flow. This statistical approach is achieved through the 
Reynolds averaging procedure which begins by decomposing the instantaneous quan- 
tities of interest (velocity and pressure) into a mean and fluctuating component. 

q(x, t ) = Q(x ) + q'(x, t) (2.9) 

where Q(x) is the time averaged mean component, i.e. 

1 r t+T 

Q(x) = ^irr^ - jf q(x , t)dt (2.10) 

and q'(x, t) is the turbulent fluctuation. Since it is not feasible to resolve the turbulent 
fluctuations, we are instead interested in determining the effect of these fluctuations 
on the mean flow. 

After decomposing the Navier-Stokes equations (Eqs. 2. 6-2. 7) and applying the 
above time average we are left with the Reynolds Averaged Navier-Stokes equations 
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for an incompressible fluid. 


dP d 

dxi dxj 


(2.11) 

( 2 . 12 ) 


where Ui,P,Sij are the mean velocity, pressure, and strain rate, respectively, and u' 
represents the fluctuating velocity. By comparing eqs. (2. 6-2. 7) with eqs. (2.11-2.12) 
it is seen that aside from the replacement of the instantaneous variables by the mean 
values, the only difference between the instantaneous equations and the time-averaged 
equations is the appearance of a momentum flux which acts as an apparent stress. 
This additional parameter 


Tij = pu'i iij (2.13) 

is referred to as the turbulent or Reynolds stress tensor. In order to solve the time 
averaged equations of motion for the mean flow properties, we must have a means of 
computing r^ . A similar procedure applied to the energy equation produces an un- 
known turbulent heat flux vector, q tj = —pu'jh'. Since the current work concentrated 
on incompressible flows the procedure for calculating the turbulent heat flux vector 
is deferred to Refs. [10, 11]). 

The calculation of a turbulent flowfleld has now become dependent on the accu- 
rate calculation of the Reynolds stress tensor. It should be noted that an equation 
for the Reynolds stress, — pu'u' can be derived (see Appendix A), by taking higher 
moments of the Navier-Stokes equations. However, when these moments are taken, 
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additional unknowns are generated (third order correlations). This inability to close 
the equation set (i.e. more unknowns than equations) is known as the ’turbulence 
closure’ problem and is due to the non-linearity of the Navier-Stokes equations. The 
purpose of turbulence modeling is to devise approximations for the unknowns in terms 
of the known flow properties in such a way that the number of equations is equal to 
the number of unknowns. 

2.2 Eddy Viscosity 

By examining the molecular transport of momentum 2 , the kinetic theory of gases 
provides an expression which relates the viscous stress tensor to the strain rate tensor 
through a coefficient of viscosity, /i. 

Uj = 2 uSij (2.14) 

= ^pVthlmfp (2-15) 

where v t h is the thermal velocity and l m f p is the mean free path. By making an analogy 
between the behavior of turbulent fluctuations and random molecular fluctuations, 
an expression for the turbulent viscous stress tensor (Reynolds stress) in terms of an 
apparent or eddy viscosity, /x t is obtained. 


= 2 HtSij (i 7^ j) 

(2.16) 

_ 1 , 
ftt — 2 P® char** char 

(2.17) 


where v char and l c h aT are some characteristic velocity and length scales representative 
of the turbulence. The problem of turbulence modeling has now been reduced to one 
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of calculating Eq.(2.17). Keep in mind that although Eqs. (2.15) and (2.17) appear 
very similar, they do have one major difference. Namely, the thermal velocity and 
mean free path are properties of the fluid, whereas the characteristic velocity and 
length scales are flowfield properties. This difference is the reason why it is possible 
to give an a priori description of p but not p t . 

2.2.1 Levels of Modeling 

There are several levels of turbulence modeling to consider ranging from the sim- 
plest algebraic model to very complex Reynolds stress formulations. A brief descrip- 
tion of each method will be included for completeness, however, readers interested in 
a more in depth analysis should see Ref. [2]. 

The first, and simplest method considered, is an algebraic turbulence model. Typ- 
ically, the characteristic velocity is given by a simple algebraic relation and the mixing 
length is based on some characteristic length of the flow. For wall bounded flows l char 
is often based on normal distance from the wall, and for free shear flows it is typically 
proportional to the width of the shear layer. Herein lies the problem for algebraic 
models. This proportionality constant varies for each of the different types of free 
shear flows. This type of model, which requires a priori information about the flow- 
field (other than initial and boundary conditions) is referred to as an "incomplete” 
model. 

The next level of modeling is termed a one-equation model. Instead of using an 
algebraic relation, an equation is solved (typically the turbulent kinetic energy) for the 



12 


characteristic velocity. This formulation still requires the specification of a turbulent 
length scale, and therefore, by analogy with the algebraic method, the one-equation 
model is also an incomplete model. 

The next level of complexity for turbulence modeling involves solving two equa- 
tions for the two characteristic scales. A two-equation formulation is desirable because 
it is the simplest form of a 'complete’ model, where complete refers to the fact that 
only initial and boundary conditions for the two-equations must be specified. No 
a priori information about the flowfield is required. Higher order ’complete’ turbu- 
lence models are available such as the Reynolds Stress model, however, solving this 
imposes a great deal of complexity since five new equations must be solved for the 
five independent components of the Reynolds stress tensor. Another reason for the 
choice of a two-equation model over a stress model is due to the fact that most of the 
problems with existing turbulence closure models have been traced to an inadequate 
dissipation rate equation 12 . When this is taken into consideration along with the fact 
that solving a Reynolds stress model requires a dissipation rate equation, the choice 
of a two-equation model becomes obvious. It should be noted, however, that once an 
improved dissipation equation is developed, an extension of the model to a Reynolds 
stress formulation becomes the obvious next step. 

Prom Eq. (2.17) it is seen that you can solve one equation for a velocity scale and 
one equation for a length scale, or alternatively, you can solve for a velocity and a time 
scale. This is explained by noting that on dimensional grounds alone, a characteristic 
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velocity and length can be combined to give a characteristic time scale. For the 
current work we have chosen to solve one equation for the turbulent kinetic energy 
(characteristic velocity scale) and one equation for the dissipation rate of turbulent 
kinetic energy (characteristic time scale). Since the TKE equation is assumed to 
be a well defined/modeled equation, and since most of the debate over two-equation 
models focuses on the dissipation equation, the current formulation uses the standard 
k equation (see Appendix B) and a new dissipation equation based on the fluctuating 
vorticity. It should be pointed out that a variety of other dissipation rate equations 
could have been considered (i.e. u, e), however, the reason an equation based on the 
fluctuating vorticity was chosen can be seen by noting some predominant processes 
occurring in a turbulent flowfield. 

Turbulent flows are dominated by three dimensional, highly rotational processes. 
Furthermore, vortex stretching is noted to be the most effective means of extracting 
energy from the mean flow and transferring it from the large eddies to the smallest 
eddies 1 . The smallest eddies, in turn, dissipate the energy into heat through molec- 
ular viscosity. This predominance of fluctuating vorticity in maintaining a turbulent 
flow suggests that a natural route in our search for a second equation (i.e. a charac- 
teristic time scale) would be through an examination of the vorticity equation. More 
specifically, it will be shown that the fluctuating vorticity can be used to calculate the 
dissipation rate as required to close our equation set (see Appendix B). A derivation 
of the Enstrophy equation is provided in Appendix D. It should be noted that the 
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Enstrophv equation has also been the subject of detailed investigations by Bernard 
and Berger 13 and Raul and Bernard 14 . It should also be noted that the current 
formulation for a dissipation equation can be easily incorporated into existing CFD 
codes by substituting the enstrophy equation for the V or ‘u>’ equations. 

As previously mentioned, the current method is based on solving one equation for 
the turbulent kinetic energy (per unit mass) 

k = iu'u' = ^ ( u' 2 + v' 2 + w 12 ') (2.18) 

and one equation for the enstrophy (vorticity variance) 

C = (2.19) 


where u\ is the fluctuating velocity and u o[ is the fluctuating vorticity. We now define 
the eddy viscosity as 


Pt 




(2.20) 


where is the structural factor (see Appendix G), p is the density, k is the turbulent 
kinetic energy, v is the kinematic molecular viscosity, and £ is the enstrophy. In order 
to determine the turbulent stresses we apply the Bousinesq approximation for an 
incompressible medium. 


2 

Tij = 2 ntSij - -pkdij ( 2 . 21 ) 

where dij is the Kronecker delta. Note that (—^pkSij) has been added to the stress 
relation (Eq. 2.16) to ensure that 


At = -pu'iii'i 


(2.22) 



15 


= —2 pk 


(2.23) 


as defined in Eq. 2.18. 

From the above two-equation formulation it is seen that only the initial and bound- 
ary conditions for k and C must be specified, therefore, this type of turbulence model 
is a ’complete’ model. 



3 Modeling Unknown 
Correlations 
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The derivation of the exact incompressible k-( equations are given in Appendices B 
and D. The final forms have been reproduced here for convenience. 


Dk 

Dt 

DC 

Dt 


Tn dU z du'i du'i d [ dk 1 p'u' 

p dxj dxk dxk dxj [ dxj 2 3 p 

-2u'u4^^ — + 2 u'iUjSjj + 2 + 2Q j uj , 1 s' ij 

d 2 Q dco 1 du)[ 

4 -v _ 2v — - — - 

dx'^ dxj dxj 


(3.1) 


(3.2) 


Now that we have our equation set the next step is to model the unknown quanti- 
ties in terms of known flow variables. The manner in which this is done represents the 
primary difference between the k-C, model and other two-equation models. Tradition- 
ally, the high turbulent Reynolds number hypothesis is employed, whereby, the largest 
scales are unaffected by the fluid’s viscosity and the smallest scales are isotropic, or 
independent of the larger scales. If the smallest scales are independent of the mean 
flow, this implies the following terms (which are dependent on the mean flow) from 
Eq. (3.2) can be neglected 


2it'u;'^ 1 , 2u' i u>' j Sij 


2 


The above terms represent the fine scale correlations which describe the detailed 
mechanism of the dissipation process and thereby control the near wall behavior. 
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Thus, these terms must be retained if one hopes to accurately predict the skin friction 
and heat transfer. Therefore, the current model retains and models each of the above 
described terms. 

In addition to retaining all of the fine scale terms in the exact equations a set of 
modeling guidelines were developed in order to impose a consistency to the modeling 
process and also in the hopes of more accurately capturing the underlying physics. It 
should be noted that these guidelines are necessary, but not sufficient conditions, for 
modeling the unknown correlations. 

1) Dimensional consistency was enforced. 

2) All modeling maintained coordinate system independence. 

3) Gallilean Invariance was imposed. 

4) All modeling was based on linear relationships ; this is consistent with Eq. 2.21 

During the modeling we encountered a number of second order symmetric and 

skew-symmetric tensors which required modeling. In order to be mathematically 
consistent the unknown symmetric tensors are modeled using known symmetric ten- 
sors, namely, r^, S tJ , and 5 tJ . Similarly, the skew-symmetric tensors are modeled 
using Qij, tij k , and ~ '§£"§£'')■ The modeling for each of the six unknown 

terms is shown in Appendix E, with the resulting equation set shown below. 
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+ ( a 3P(bi-j + -Sijp£ ) Sij — 
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2 S^lqsuci + !hg.n i n j s ij + 

kv U 2 


kn 23gi " m ( t ) 


9/c d£ \ Qj 


<9x( <9x r 


i 2 


s 2 

(3.5) 


1 + 


<7 


a 2 At. 
P~ ^S 2 


[£ / at/, \ 
2 W*J 


where 


i~ij 2/i(iS < ij ^$ijpk , 

S 2 = S tJ S tJ , 0 2 =(1A 


k 


/2fc = 


J'v/C 




2 

it 


Now that we have the modeled equation set, all that is remaining is to deter- 
mine the appropriate boundary conditions and calculate the unknown closure coeffi- 
cients (q 3, /?4j /?5j /?6, /?7 > /^8> &Q j <2pi * 
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4 Boundary Conditions 

With the current formulation for solving turbulence we have introduced two ad- 
ditional unknowns to our equation set, namely, k, and £. A wall and free stream 
(subscripts w and oo, respectively) boundary condition is required for each of the 
added variables. 

At the wall it can be shown that k and its derivative are zero and C is finite (see 
Appendix I). Therefore, we have 


k w = 0 


dk 

drj 


= 0 


(4.1) 

(4.2) 


C has no physical boundary condition at a solid surface. 

For the current work we have examined two boundary conditions for £ w . The first 


was a simple first order extrapolation, i.e. we set 


<K 

dr} 


W 


= 0 


(4.3) 


where r} is the normal coordinate from the surface. The second boundary condition for 
was developed from Eq. (4.2) and the k equation. In the near wall region both the 
convection and production terms are negligible, and therefore, for an incompressible 
flow Eq. (3.1) reduces to 


d_ [1 dk 
dy 3 dy 


(4.4) 
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If we define a one-sided difference as 


dq _ q i+ 1 ~ 
dy Vj + 1 - 


(4.5) 


and we use a subscript convention of 1 and 2 to denote the the first and second points 
off of the wall, respectively, then we end up with 


1 

dk 

dk 


3 Ayx 

dy 

x dy 

U)- 


(4.6) 


or, using Eq. (4.2) 


1 dk 
3Ayi dy l 

(fc 2 — Ag) 
3AyiAy 2 


(4.7) 


where, 


Ar/x =y 2 ~yi 


(4.8) 


The merits of both boundary conditions will be discussed in the results section where 
the simple extrapolation is labeled as bc\ and the above formula (Eq. 4.7) is labeled 
as &c 2 . 

Now, all that is remaining is to specify the free stream values for each of the three 
variables. This is done by specifying a turbulent intensity, T and a ratio of turbulent 
to laminar stresses, and solving for the remaining variable, Coo- 

The turbulent intensity is defined as 


r = 


K[ 

u, oo 


(4.9) 
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and if you assume the flow is isotropic ( u ' = v' = w') in the free stream k can be 
written as 


kno — ~ a 


a 


(4.10) 


Eqs. 4.9-4.10 can be combined to give 


koo = 2 (TUooY 


(4.11) 


The turbulent intensity is often provided by experiment and is typically O (1%). 

It should be noted that as you leave the boundary layer and approach the free 
stream all variables, k , £, and decay to small values. Transients during the solution 
procedure can cause large spikes in the ^ ratio since fi t is defined as a ratio of In 
order to prevent numerical difficulties resulting from dividing by a very small number, 
a bound is placed on the allowable values of This limit is derived by examining 
the isotropic decay of homogeneous turbulence (see Appendix K) and is given as 


Lh 


<C "((^2))' 


where C^, /?s, and 6 are constants (see Table 9.3) 



5 Evaluation of Closure 
Coefficients 


The final step in the development of the A:-£ model is the determination of the 
unknown closure coefficients. This procedure was accomplished in four stages: a log- 
law analysis, an examination of similarity solutions for free shear flows, a boundary 
layer solution of a fiat plate and its wake, and a number of Navier-Stokes solutions for 
more complex flowfields. A brief description of each method will be provided along 
with an explanation of why the method was chosen. 

5.1 Log- Law Analysis 

The log-law refers to a region of the boundary layer over which a logarithmic 
velocity profile is valid. If we examine the /:-£ equation set in this region, and ensure 
that the correct scaling is met, the following relation is obtained. 


<7 C = 


6/c 2 


(5.1) 


[3/^ - 3 a 3y /C^ - 4A - 8ft] 

This has allowed us to express one of our unknown closure coefficients in terms 
of the others, thereby, reducing the amount of unknowns by one. A more in-depth 
explanation of the log-law region, as well as a derivation of Equation (5.1) is provided 
in Appendix G. 
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5.2 Similarity Solutions 

A number of model constants were chosen in order to optimize the solutions for 
a variety of free shear flows. The free shear flows considered for the current research 
include wakes, mixing-layers and a variety of jets (planar, round, and radial). 

This method of determining model constants was chosen for several reasons, the 
first of which being, free shear flows, by definition, are not bounded by solid surfaces, 
and therefore, the complications imposed by solid boundaries can be avoided during 
the initial evaluation of the turbulence model. These complications include, amongst 
others, a vanishing turbulent Reynolds number and the lack of physical boundary 
conditions for the dissipation equation (£) at a solid surface. 

Another reason free shear flows were chosen was due to the fact that, for each 
of the flows under consideration, a self-similar (or self-preserving) mean flow state 
is attained far downstream of the generating body. This self-similarity implies that, 
at a sufficient distance downstream, the profiles become invariant if plotted against 
a similarity variable, rj. In other words, the equations are transformed from an x- 
y to an x-r] coordinate system (where 77 = yf(x)). Then, for the types of flows 
under consideration, it is possible to eliminate the x-dependence, thereby reducing the 
partial differential equations to ordinary differential equations (ODE). These ODE’s 
are not computationally intensive, and therefore, a large number of model constants 
can be evaluated with very little effort. 

An additional reason for choosing free shear flows to calibrate the turbulence 
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model is because of the vast database of experimental data. This abundance of 
turbulent data for free shear flows allows one to easily evaluate a given model for 
a wide variety of flowfieids. It should be noted, however, that there was a lack 
of consistency amongst the various experiments in measurements of the turbulent 
kinetic energy and the normal stresses. Therefore, the constants were optimized by 
emphasizing a comparison with measured growth rates and maximum shear stresses, 
which were in close agreement (within 10%) in the various experiments. 

The final reason free shear flows were used as a starting point for determining 
closure coefficients stems from the consideration that if the new model is incapable of 
predicting growth rates and shear distributions of all free shear layers, then it can not 
have any advantage over existing k-e and k-u> models and, as such, merits no further 
development. 

It seems logical that if a model can not predict simple free shear flows, then it 
would have difficulty in accurately predicting a turbulent flowfield over a complex 
vehicle. Therefore, the obvious place to begin an optimization of a generalized tur- 
bulence model was through the examination of a variety of free shear flows. 

5.3 Boundary Layer Solutions 

The next step in determining the model constants is to test and refine the model 
for a wall bounded flow. A boundary layer solution 15 was chosen for this task due 
to its computational simplicity. This is a result of the fact that boundary layer code 
run times are much smaller than for a Navier-Stokes solver. A brief description of the 
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methodology used in the boundary layer code is provided in Appendix H. Because 
boundary layer codes cannot handle separated flows, further development required 
the use of a Navier-Stokes solver. 

5.4 Navier-Stokes Solutions 

Up to this point we have considered a variety of methods which make a number of 
simplifying assumptions about the flowfield in order to reduce the complexity of the 
Navier-Stokes equations. However, there are many flowfield features which cannot be 
examined using the simplified methods discussed up to this point. These include, but 
are not limited to, stagnation points, shockwaves, large adverse pressure gradients, 
and separated flows. In order to assess the models performance in and around these 
types of flow features, a Navier-Stokes solver must be employed. A description of the 
two Navier-Stokes solvers 10, 16 employed in the current research will be included in 
the results section. 

Initially, the Navier-Stokes solver was used to examine a variety of two-dimensional 
airfoils, the purpose of which was twofold. First was the development and evaluation 
of a pressure gradient term, and second was an assessment of the A;-£’s numerical 
stability. Once the development and validation was completed for the two-dimensional 
cases, the performance of the model for a three-dimensional supersonic Cylinder-Flare 


was examined. 
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6 k-( Modifications 

The previous section described how the closure coefficients were evaluated, how- 
ever, it is important to recognize that this was an iterative procedure. As the research 
progressed from the simpler to the more complex flows it was often discovered that 
a particular choice of modeling was inappropriate. Terms were re-modeled and the 
procedure was repeated until we reached the current state of the model. The modi- 
fications of the original model warrant special attention, and therefore, a discussion 
of each has been included below. The final form of the modeled k-C, equation set, 
including all of the modifications discussed below are shown in Appendix L 

6.1 /?5 Term 


The /?5 term is one of the dominant dissipation terms in the enstrophy equation. 
Originally it was modeled as 


fl»C* 

Rk 


( 6 . 1 ) 


for the free shear layers. However, when examining wall bounded flows R * is zero at 
the wall. It should be noted that for the current work the finite volume approach 
was used so this term is never actually evaluated at the wall. However, as the wall is 
approached Rk can become extremely small, thus causing numerical problems. This 
numerical difficulty is a result of allowing the turbulence time scale to be less than 
Kologomorov’s time scale. Calculation of Rk using Kolomogorov’s microscales 1 yields 
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a value of about unity, 0( 1). To enforce the requirement that the turbulent time scale 
can not be less than Kolmogorov’s time scale, the dissipation term can be written as 


fed 

Rk + S 


( 6 . 2 ) 


where (5 = 0(1). Determination of S will be explained in the next section on wall 
damping functions. 


6.2 Wall Damping Functions 


For high turbulent Reynold’s numbers the eddy viscosity is defined as 


A k = 


K 


(6.3) 


For wall bounded flows the eddy viscosity and the turbulent Reynolds number are 
both zero at the wall. Because of this the above expression is traditionally multiplied 


by a wall damping function, / M . 


IM = 


pC»k 2 U 

K 


(6.4) 


There is no unanimity on the form of / M in the near wall region 17 . Initially 18 the k-Q 
model was implemented using the following form for the wall damping function 


/„ = min (/,), 1.0) 



(6.5) 

( 6 . 6 ) 


An examination of the above expression reveals that this term is unity except for 
a very small region near the surface ( y + < 100). The problem with employing a wall 
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damping function is the introduction of geometrical dependence into the model. The 
’y’ is the normal distance to the nearest wall, however, for complex three-dimensional 
flows its calculation is somewhat cumbersome. Therefore, in order to make the k -( £ 
model more amenable to three-dimensional solutions, the wall damping function has 
been eliminated in recent implementations[19]. 

During the modeling process (see Appendix E) the following 




(6.7) 


was re-written using a symmetric and a skew-symmetric term with the symmetric 
portion being modeled as 


v_ t_ 

2a r 


<9Q, dfi, 

+ 3 


dxj dxi J \dxj 


dn 


( 6 . 8 ) 


The above term does not contribute to the similarity solution of free shear flows, 
and as a result, a T was assigned a value based on log-layer considerations which 
was inappropriate. Later, it was shown that the value of o r was somewhat depen- 
dent on the model constants that appeared in the damping function, A previous 
investigation 18 showed that it was possible to eliminate the a r term by simply modi- 
fying the constants. This result raised an interesting question: if one can develop 
an f u that eliminates Eq. (6.8), then can one eliminate / M and keep the above term. 
It turned out that this was indeed possible by choosing the model constants a r and 
6 as 


S = 0.1 , a T — 0.07 


(6.9) 
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From our previous analysis it was shown that 5 = 0( 1). A major consideration 
that entered into the selection of 5 is transition from laminar to turbulent flow. In 
the absence of a transition model, small values of 6 are required to allow transition 
to take place. 

6.3 Symmetric Portion of — 

The symmetric portion of this term is written as (See Appendix E) 


€ mij 

d (u'mU'l) 

dk j 

2 

dxi 

i 

£ 

3 

1 


An evaluation of Eq. 6.10 has shown its impact to be negligible for all of the cases 
considered. Furthermore, an examination of this term in a three-dimensional coor- 
dinate system reveals that you are required to calculate the derivative of each of 
the Reynolds stresses in each direction. This can be quite costly computationally, 
and therefore, in order to simplify the k-^ model this term has been neglected. No 
noticeable change in results was noticed from the omission of this term. 

6.4 Pressure Gradient 

In eddy viscosity models such as the k-e and k-uj, the production of turbulent 
kinetic energy is related to the mean strain rate. This is large in regions of significant 
curvature, and therefore, eddy viscosity models typically generate excessive produc- 
tion of k in these regions. This in turn increases the value of which delays or 


inhibits separation. 
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The successful models that were developed to address separated flows were aimed 
at mimicking conditions that exist in separated regions by reducing eddy viscosity or 
turbulent energy production. Thus, the Johnson and King 8 model was designed to 
reduce the eddy viscosity in the wake region of boundary layer flows characterized 
by adverse pressure gradients. Reduction of the eddy viscosity is also the motivation 
behind the shear stress transport model of Mentor 20,21 . For the current model a 
reduction of eddy viscosity in adverse pressure gradient regions can be imposed in 
one of two ways. First, one can model the term 

1^1 (611) 
that appears in the k equation (Eq. 3.1) in terms of a mean pressure gradient. This was 
attempted by Rao and Hassan 3 but proved to be unsatisfactory. Alternatively, one 
can develop a term in the enstrophy equation based on the mean pressure gradient. 
Although the pressure does not appear explicitly in the exact (incompressible) £ 
equation, the effect of the pressure is noted by examining the term This term 
can be expressed as a second derivative of mean velocity which is related to the mean 
pressure gradient through the momentum equation. Thus, the £ equation can have 
an explicit dependence on the mean pressure gradient. Based on this, the following 
term was originally added to the k-£ model 



This term produces £ , and therefore, reduces the eddy viscosity in an adverse 
pressure gradient flow as desired. However, the problem with a term of the above 
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form is that its always positive, and therefore, always dissipative of p t , regardless 
of whether favorable or adverse pressure gradients are considered. The result was 
that the above term was excessively dissipative around the leading edge of airfoils. 
In order to correct this problem we were looking to develop a term which took into 
account the sign of the pressure gradient (i.e. favorable or adverse). This was done 
by re-modeling the above term as 

!_ (6 . 13) 

dx, i/P(Tp (1 + <5p) 1 ' 


where 


kR t ( dp \ 

C \dxij 



(6.14) 


This term (Eq. 6.13) represents the derivative of pressure along a streamline, and 
therefore, is negative (produces p t ), and positive (dissipates p t ) for favorable and ad- 
verse pressure gradients, respectively. The (1 + Sp ) is a compressibility correction and 
has been added to prevent the pressure gradient term from becoming excessive around 
strong shocks. Note that, in the absence of mean pressure gradients, the modeling 
reverts to that appropriate for constant pressure solutions. Thus, no adjustments 
in the previously calibrated model constants are necessary. It should also be noted 
that there are a number of possible methods for modeling the pressure gradient term. 
Therefore, Eq. 6.13 should be considered a tentative term until the efficacy of other 


forms can be examined. 



6.5 Compressibility and Transition Modifications 


Modifications to the k-£ model which take into account compressibility effects are 
currently being developed by Alexopoulos 11 ’ 22 . He has developed a number of com- 
pressibility terms (see C l5 terms in Appendix L), all of which have been included 
in the current results. These additional terms have had very little effect, as expected, 
since the cases currently being considered were all either subsonic or transonic, or as- 
sumed adiabatic walls (Cylinder-Flare). Also, a transition model has been developed 
by Warren 23, 24, 25 based on the k-£ model. This transition model was not used for 
the current results. For an in depth discussion of both compressibility and transition 
and how they relate to the k-( model the reader is referred to Refs. [11] and [25]. 

6.6 Numerical Stability Modifications 

The current model, like all turbulence closure models, requires some modification 
to improve its numerical stability. The first stability improvement is made by setting 
a minimum £ in the Navier-Stokes codes. The reason for this minimum can be seen 
by examining the typical behavior of k and £. The enstrophy is finite at the wall, 
attains a maximum in the boundary layer, and decays to some small value outside 
of the boundary layer. Similarly, k is zero at the wall, reaches a maximum in the 
boundary layer, and then decays to some small value outside the boundary layer. An 
examination of /x t 
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shows that since both k and ( are rapidly approaching very small values at the edge of 
the boundary layer, numerical instabilities can arise. In order to improve the stability 
and prevent oscillations in /j, t a floor is typically set for the dissipation, £. For the 
current work a minimum £ is set as 

C = max (c, (6.16) 

Minimums as large as ^ were evaluated and shown to have no effect on the final 
converged solution. 

Two additional modifications were made in order to improve the numerical sta- 
bility of the model. Namely, both the pressure gradient and /3 8 terms are maxed with 
zero (See Appendix L. This forces both terms to be shut off in regions where they 
dissipate £ (produce For the pressure gradient term this was done because this 
term produced excessive n t around the leading edge of the airfoils. It should be noted 
that this is only a transient phenomena which occurs during the convergence of the 
solution. Once the solution converges the max can be removed without adversely 
effecting the final solution. A similar procedure applies for the /? 8 term. 

6.7 Summary of Modifications 

It should be noted that as modifications to the model are made one underlying 
assumption must hold. Namely, that the modification does not affect the previously 
considered cases. In other words, the addition of the 6 to the fa term to improve the 
results for the wall bounded flows can not affect the free shear results. Similarly, the 
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addition/removal of a wall damping function can not affect either the boundary layer 
solution of the flat plate or the free shear results. This assumption holds for each of 
the additional modifications made in this section. Also, as additional modification are 
made which take into account more complex flowfield features (e.g. compressibility, 
combustion, etc), this assumption must remain valid in order to develop a generalized 
turbulence closure model capable of solving complex three-dimensional flows. 

The final form of the A>£ model with all of the aforementioned modifications, 
including the compressibility terms developed by Alexopoulos 11, 22 , can be seen in 
Appendix L. 
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7 Results 

7.1 Free Shear Layer 

For the free shear layers we have compared the current model to both experimental 
data and the results from the k-e turbulence model. The results from the k-ui model 
were not included due to the fact that the free stream boundary condition for u> must 
be adjusted for each of the different types of shear flows. 

We begin with the results for the self-similar incompressible planar wake. A 
schematic of this flowfield is shown in Figure (10.1). Figure (10.2) shows the defect 
velocity, W (non-dimensionalized by centerline defect velocity W 0 ) versus y (non- 
dimensionalized by wake width). For this case we have compared with a variety of 
profiles at various | locations from the experiments of Patel 26 as well as with the k-e 
turbulence model. Note that the experimental data shows that a self-similar velocity 
profile is obtained at an | of approximately 79. The k-C, model shows excellent 
agreement with experiment and also shows a noted improvement over the k-e model. 
Figure (10.3) shows a similar comparison for the shearing stress, r along with the 
experimental data of Pot 6 . Note, that two sets of data were included at an | of 960 
in order to show the level of asymmetry in the experimental data. Also included was 
a comparison with the asymptotic solution obtained by assuming a constant eddy 
viscosity, fi t . Once again the k-£ model shows excellent agreement with experiment 
as well as with the asymptotic solution. Figure (10.4) shows the turbulent kinetic 
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energy profile, and while there was no available experimental data for this case, it 
has been included for completeness. The effect of the improvements in the velocity 
and stress profiles can also be seen in Table 9.1 which shows a percent error in the 
growth rate of 14% and 30% for the k-C j and k-e models, respectively. 

The next set of Figures (10.6-10.8) show the comparisons for velocity, shear stress, 
and turbulent kinetic energy for the two-dimensional planar jet (see schematic in Fig- 
ure 10.5) compared with experiment 27, 28, 29 and the k-e model. The results between 
the k-( and k-e models are essentially identical except for the formers slight over pre- 
diction of shear stress. Both models agree extremely well with experimental profiles 
and, therefore, both do a very good job of predicting the growth rate (Table 9.1). 
Bradbury 28 indicates that his measured spreading rate is not exactly proportional to 
x. The departure from true self-preservation, however, is so small than no adjustment 
is necessary in Fig. (10.6). 

Figures (10.9-10.11) show the comparison of mean velocity, shear stress, and tur- 
bulent kinetic energy for the round (antisymmetric) jet along with the results of the 
k-e model and experiment 30, 31 . The k-C model does an excellent job of predicting 
the velocity and shear stress profiles. Note, however, that the centerline TKE is un- 
der predicted. This is explained by noting that during the optimization process we 
were primarily concerned with accurately reproducing growth rates, velocity profiles, 
and peak shearing stresses. The primary reason for considering the TKE prediction 
to be of secondary importance can be seen by noting that there was a large degree 
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of variation in the measured TKE for the various experiments, and therefore, these 
measurements were deemed to be less reliable. Therefore, we primarily concentrated 
on matching velocity and shear profiles as well as growth rates. Note that the k-Q 
models prediction of growth rate lies within the experimental range while the k-e 
shows a 26% error if compared with the upper bound (Table 9.1). 

Figures (10.12-10.14) shows a similar comparison for the radial jet. Although no 
experimental data was available, other than growth rates, the plots were included for 
completeness and to show the similarity between the k-C > and k-t models. Available 
data 32, 33, 34 suggests that profiles for the radial jet approach those for the plane jet. 
Since both models were in good agreement with experiment for the plane jet, the 
indicated agreement in Figures (10.12-10.14) is expected. Note that both models did 
an excellent job of predicting the growth rate for this case. 

The final set of figures (10.16-10.18) show the comparisons for the mixing layer 
(see schematic in figure 10.15). In these figures 77 is defined as 


V = - 


(j)-(iL 
(D„, - (Do, 


(7.1) 


where , etc., denotes the locations where the ratio of the mean velocity to the 

free stream velocity is 0.1, etc. For this case the mean velocity is accurately predicted, 
however, the shear stress and TKE both failed to reproduce the experimental data. 
The growth rates were within 5% for the k-£ model compared to 15% for the k — e 


model. It should be noted that the theory shows good agreement with the predicted 
shear stress in the streaming side of the measurements but not with the zero velocity 
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side. Patel 35 points out that measurements with a normal hot-wire probe became 
unreliable in the zero velocity side. However, no estimate of inaccuracies were given. 
Thus, a more accurate measurement in the zero velocity side is needed before a 
definitive statement can be made regarding predictions of the theory in this region. 

A summary of the calculated spreading rates is shown in Table 9.1. The spreading 
rates are defined in a manner consistent with that given in Ref. [2], Moreover, the 
e comparison values were obtained from Table 4.2, Ref. [2]. In general, computed 
growth rates represent a noted improvement over those obtained from the k — e or 
k — cj models. 

After optimizing the free shear flows and obtaining a tentative set of closure 
coefficients the next step was to examine the turbulence model for a wall bounded 
flow. This was done by examining the boundary layer solution of a flat plate and its 
wake. 

7.2 Flat Plate 

The next step is to calibrate the k-Q turbulence model for a wall bounded flow. 
As was previously discussed, the original form of the model contained a wall damping 
function with two unknown closure coefficients (C^,, C^). In order to evaluate these 
constants an incompressible flow (M = 0.01) over a flat plate was considered. 

Figures (10.19-10.21) show plots of k + , £ + , and vs. y + for the near wall region. 
The current model is compared with the experimental measurements of Schubauer 4 
and Laufer 5 as well as the ’’average” of the available data 17 . As is seen from the figures, 
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the results compare well with the results summarized in Ref. [17]. From Figure (10.19), 
it is seen that the peak value of k + has been slightly under predicted, however, the 
current results occur within the range typical of two-equation models (k + « 2.8 — 4.0). 
Figure (10.20) shows that both £ + and e + have a ’’peak” near the wall. It should be 
mentioned that the e + is not from a solution of the k-e model, rather it is determined 
by comparing the k equation for the £ model to that of the standard e model and 
noting that 


e = K + 


2v d 2 k 
3 dydy 


(7.2) 


Therefore, the e + shown, is that calculated by using the above expression in conjunc- 
tion with the solution of the k-( model. The current model accurately reproduces the 
experimental data of Laufer, however, there has been some recent controversy over 
the validity of this data in the very near wall region. The DNS data for this case 2 
shows no such ’’peak” in the near wall region. 

The constant, ’B’ from the log-law relation 


u + = — In (y + ) + B 
k v ' 

is typically taken to be 5.0. A plot of u + — Hny+ shows that ’B’ is not exactly 
constant, but rather, varies slightly from approximately 4.4 to 5.2 over a y + range 
from 50 to 1000. 

The sublayer behavior of the TKE is shown in figure 10.23. Since n (k = k a y n ) 
could not be evaluated at y = 0, a polynomial curve fit was used to obtain a value of 
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2.08 at y = 0. This is in very close agreement with the theoretical value of 2.0 (see 
Appendix I). 

Figure (10.24) shows the skin friction coefficient (based on edge values) vs. the 
vanDriest 36 correlation. Excellent agreement with the current model is shown. 

Figure (10.25) shows the wall damping function. Note that this parameter is 
non-unity only in a small region very close to the wall ( y + < 80). 

Once acceptable results were obtained for the flat plate, the boundary layer solu- 
tion was used to examine the near wake properties and as a means of validating the 
self-similar wake results. This was done by marching off the end of a flat plate and 
continuing the solution procedure downstream until similarity profiles were attained. 

Figs. (10.26-10.30) demonstrate the fc-^’s prediction of a variety of near wake 
parameters including growth rate, centerline defect velocity, and peak shearing stress. 
The results of the current model are compared with the experimental data of Pot 6 
and Weygandt and Mehta 7 . Excellent agreement is indicated. Also shown are the 
stress and velocity profiles in the far downstream region. These figures show that a 
self-preserving state was attained and that the marching procedure has reproduced 
the similarity results within a few percent. 

It should be mentioned that all of the figures considered up to this point used 
the wall damping function, f fl . It should also be noted, however, that elimination of 
the wall damping function made no noticeable changes in the flat plate results. This 
was the criterion used for the elimination of the damping function. In other words, 
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the model had to give the same results with or without the / M . All of the results 
for the flat plate without the damping function were within a couple of percent, and 
therefore, they have not been reproduced here. Also, note that the elimination of 
the fp has no effect on the free shear flows since was unity for all of those cases. 
All of the remaining results were obtained after the elimination of the wall damping 
function. 

7.3 Homogeneous Shear Flows 

This section examines the predicted decay rates of turbulence in a homogeneous 
shear flow. It should be noted that even though the current model constants were 
not chosen by considering the decay of homogeneous turbulence, close agreement 
with the k — t model, whose constants were chosen to give the correct decay rate for 
homogeneous turbulence, is shown in figures 10.32-10.33. This test case was examined 
per the suggestion of Reference [37] in order to determine if the current model would 
predict an acceptable value of the homogeneous shear parameter and also to ensure 
that the h-C, model would predict equilibrium (or constant steady state) values for 
the shear parameter and for the Reynolds stress anisotropy, 6 12 . Measurements and 
DNS 38 suggest that the shear parameter reaches an equilibrium value between 5 and 
6. Figures (10.31-10.32) show that the current model does in fact reach a constant 
steady state value and also does predict an acceptible homogeneous shear parameter. 

Figures (10.33-10.34) show the comparison of the h-C, and the k-e models in the 
prediction of k + and c + . These figures represent the time decay of the turbulent 
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quantities and show very similar results between the two models. One final note 
should be made regarding homogeneous shear flows. Namely, that for this case the 
fc-C model gives 6n = b 2 2 = 633 = 0, which has been shown to be incorrect [39]. 
This is not a limitation of the current two-equation model, rather, it is a limitation 
in the current Bousinesq approximation. In order to be able to accurately predict 
the Reynolds stress anisotropies, the Bousinesq approximation must be modified to 
include higher order terms. 

7.4 Airfoils 

For the airfoil cases a comparison is made with the k-u model of Wilcox 2 , the 
Johnson-King 8 model, as well as with experimental data where available. The k-u 
model was chosen because it is widely known to have good agreement with experiment 
for attached wall bounded flows. The J-K model was included because it is widely 
known to accurately predict two-dimensional separated airfoils. The J-K results were 
obtained from the calculations of Rumsey and Anderson 40 . Mentor 20 made an issue 
of the sensitivity of the k-u model to free stream conditions. In all of the calculations 
presented here, k <*, and were selected and corresponding values of Coo and u^ 

were computed. We believe this is the correct way to compare the two models. 

The evaluation of the A;-£ model for two-dimensional airfoils begins by first ex- 
amining the effect of the two boundary conditions. Recall that bc x is a simple of 
extrapolation of C to the wall, whereas, bc 2 imposes a k = k 0 y 2 behavior to determine 
a wall value of C- The first case considered is a symmetric NACA 0012 airfoil at a 
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relatively low Mach number. This provides a simple flow free of both shockwaves and 
separated regions. Note from figures 10.35-10.36 that both boundary conditions give 
nearly identical results except for the skin friction, C/ e in a small region around the 
trailing edge. This is explained by noting that the first boundary condition (simple 
C extrapolation) gives a near wall k variation very close to the theoretical value (see 
Appendix I) of n = 2. Thus, imposing a y 2 variation of the near wall k has little to 
no effect, as expected. However, this is not the case for the more complex flowfields. 
Figures (10.37-10.38) show the same comparisons for the NACA 0012 run with a 
transonic Mach number. For this case small variations in C p can be seen around the 
shock, however, large variations in Cf e are noted. When using simple extrapolation, 
bc\ for more complex flows, large oscillations in the skin friction coefficient can occur, 
therefore, for the more complex flows bci has proven to be inadequate around shocks 
(a similar behavior was noted near separated regions). It should also be noted that 
bc\ proved to be much less numerically stable around shocks, separated regions, and 
stagnation points. Therefore, for the remainder of the results section, only bc 2 will 
be considered. 

The next step in evaluating the k-C, model was to ensure that grid independent so- 
lutions were being obtained. A grid study was performed for an RAE 2822 transonic 
airfoil with a small separated region. Figures (10.39-10.40) show the pressure distri- 
bution and skin friction for all grids considered. Both the pressure and skin friction 
distributions have proven to be grid independent. The results are nearly identical, 
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except for a slight variation in skin friction on the 143x51 grid. From this grid analy- 
sis it has been determined that the k-( model requires a minimum of approximately 
5 points within the sublayer (y + < 10). The 321x91 grid had an initial yf of approxi- 
mately 0.2, the 231x51 used yf ~ 1.8, and the 143x51 used ^ 3.0. Note, however, 
that the deviation in skin friction results were caused not by the initial y + spacing, 
but rather, is was primarily caused by the amount of resolution of the sublayer. In 
other words, as long as you have at least 5 points in the sublayer, you may begin with 
initial yf values of 2 or 3 and still obtain accurate solutions. This is important, and 
is being discussed in detail, because of the obvious influence on rate of convergence. 
Increasing the distance of the first y + point off of the surface can drastically increase 
the rate of convergence of the solution. Note, that for all of the remaining airfoil 
cases, an initial y + between 1 and 2 and approximately 8-10 points in the sublayer 
were employed. 

Now that the boundary conditions and basic grid requirements have been deter- 
mined the next step is to examine the results for a variety of two-dimensional airfoils 
beginning with a NACA 4412. The airfoil geometry is shown in figure 10.41. This 
case is a low subsonic (M = 0.2) stalled airfoil (a = 13.87°) and was used to examine 
the k-£'s performance on an airfoil with a large separated region but in the absence 
of shockwaves. Comparisons are made with the available experimental data 41 as well 
as with the standard k-w model. Figure (10.42) shows that the pressure distribution 
is well predicted by both the current and k-u> models. Although no skin friction data 
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was available, the experiments of Coles and Wadcock 41 showed separation beginning 
between ~ of .75- .80. Figure (10.43) shows that the k-( model predicts separation at 
c = 0.81, in excellent agreement with experiment, while the k-u model predicts de- 
layed separation and a somewhat smaller separated region. Figure (10.44) compares 
the mean velocity profiles for a variety of ~ c stations throughout the separated region. 
Again, prediction of the current theory is more consistent with the experiment. 

The next set of comparisons (figures 10.45-10.47) involve a NACA 0012 (geometry 
shown in figure 10.45) at an angle of attack, a of 8.34° and a free stream Mach number, 
Mqq = 0.55 compared with the experimental data of Harris 42 . All models predict the 
pressure distribution well . An absence of skin friction data prevents a meaningful 
discussion of which model is more accurate, however, it should be noted that both 
the current and k-u models predict only one separated region while the J-K predicts 
two. Figures (10.48-10.49) show the same airfoil at an angle of attack of 2.26° and a 
transonic Mach number of 0.799. Both the current and J-K models give an accurate 
prediction of the pressure distribution and show similar results for the skin friction. 
The k-u model predicts delayed separation which is not consistent with experiment. 

Figures (10.50-10.54) show the results for the RAE 2822 (geometry shown in 
figure 10.50). For cases 9 and 10, we have compared with the experimental data of 
Cook et. al. 43 ; no definitive angle of attack or Mach number corrections to account for 
wall interference were given. As a result, the flow conditions assumed for these cases 
were those used by Rumsey and Vatsa 44 . For case 9 (M^ = 0.73, a = 2.80°, R e = 
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6.5il0 6 ), the flow is attached and the k-u model does a better job of predicting the 
shock position and the pressure distribution. Figure (10.52) shows that all models 
do a fair job of predicting the skin friction for this case. The next case (10) involves 
shock induced separation. Both the A>£ and J-K do a good job of predicting both the 
shock position and the skin friction. Note that the k-uj model predicts both a delayed 
shock and a milder separated region. 

7.5 Cylinder-Flare 

The Cylinder-Flare case was considered in order to examine the k-C, models ability 
to accurately predict a complex three-dimensional supersonic flowfield. The geometry 
for this case is shown in Fig. (10.55) along with the test conditions. A schematic of 
the flare, along with the appropriate dimensions, are shown in Fig. (10.56). Also, in 
order for the reader to get a general idea of the flowfield being considered, a sketch 
of the postulated flowfield was provided in Fig. (10.57). This figure was taken from 
Ref. [9]. From this figure it it seen that at the 0 — 0° plane, there is a small separated 
region between the intersection of the cylinder and the offset flare. This separation 
bubble turns the flow prior to reaching the juncture, and therefore, causes the oblique 
shockwave to move upstream of the flare. Also, as you consider planes with increasing 
0, the extent of the separation region grows, and therefore, the amount of upstream 
influence is increased. This is due to the fact that as you progress from 0 = 0° to 
0 = 90° the flow becomes increasingly three-dimensional in nature due to the crossflow 


caused by the offset flare. 
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A few notes about the solution procedure are warranted. First, in order to min- 
imize the computational time required for this case the problem was split into two 
sections: an axisymmetric cylinder and a three-dimensional flare. A semi-infinite 
cylinder was computed assuming axisymmetric flow. At a given distance along the 
cylinder the jj-, C/, and 5 parameters were found to match the experimental results. 
The profiles for velocity, temperature, pressure, k, and C from this point on the cylin- 
der were then used as initial profiles for the three-dimensional flare. Note, that for 
the full three-dimensional case, a small cylinder (10 cm length) was run upstream of 
the flare so that the initial profiles would not be affected by the upstream influence 
generated by the flare. In other words, instead of running the full cylinder (90 cm) 
in the three-dimensional code, only a smaller portion (10 cm) was considered. The 
remaining 80cm was run using an axisymmetric assumption to reduce computational 
time. The above procedure provided for a significant reduction in the number of grid 
points required. 

The solution procedure used was a four stage Runge-Kutta, explicit, central differ- 
ence (standard Jameson Damping) code developed by Baurle 16 . The grid contained 
87 x 99 x 99 points in the x,y,z directions, respectively and was run with an initial 
Hi of 0.3 and approximately twelve points in the sublayer ( y + < 10). Comparisons 
of the current model are made with the k-e model as well as the data of Wideman et 
al 9 . The k-e results were obtained from Ref. [45]. 

Figures (10.58-10.60) show the pressure ratio for the & = 0°, 90°, 180° planes, 
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respectively. The pressure ratio is equal to the surface pressure over the freestream 
stagnation pressure. For the 9 = 0 plane both models accurately predict the initial 
pressure rise upstream of the juncture. This indicates that, for this plane, both 
models are accurately predicting the extent of the separation bubble. Downstream of 
the juncture the k-C models shows a pressure rise which is in closer agreement with the 
experimental data. Both models predict the asymptotic pressure region accurately. 
For the 9 = 90° plane both models again do a good job of predicting the initial 
pressure rise. Aft of the juncture both models overpredict the pressure ratio, but they 
both return to acceptable levels in the downstream region. The 9 = 180° plane once 
again shows that the k-C model does a slightly better job of predicting the pressure 
ratio upstream of the juncture. Aft of the juncture the k-C models shows a slightly 
larger peak, however, it returns to an acceptable level in the downstream region. 
From the three previous figures it should be noted that the upstream influence is a 
minimum at the 0 = 0° plane where very little crossflow velocity exists. It continually 
increases until the 9 = 90° plane is reached at which point the upstream influence 
is a maximum (noted by the distance between the juncture and the initial pressure 
rise). This value then remains relatively constant until the bottom symmetry plane 
is reached at 9 = 180°. 

Figures (10.61-10.63) show the skin friction coefficient (non-dimensionalization 
uses free stream quantities) for the 9 = 0°, 90°, 180° planes, respectively. For the 
9 = 0° plane the k-C, model does a better job then the k-e in predicting Cf downstream 
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of the juncture. The k-e model overpredicts the asymptotic value for skin friction. 

The Cf for the 9 = 90° plane shows that both models mimic each other very closely, 
however, they also significantly overpredict the skin friction data in the downstream 
region. The validity of the experimental data for this plane has been called into 
question by one of the co-authors of the original paper 46 , and therefore, a meaningful 
discussion of the results is not possible. It should be noted that even though the skin 
friction data for this plane is being questioned, the pressure data remains valid 46 . 

The final figures shows the skin friction comparison for the 9 — 180° plane. The 
model slightly overpredicts the separated region, however, the model does an 
excellent job downstream of the juncture. 

In summary, the A;-£ model without damping functions has performed as well as 
the k-e model which employs damping functions. Overall, good agreement is indicated 


with the current model. 
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8 Concluding Remarks 

The current A>£ model has proven to be a noteworthy step towards the develop- 
ment of a generalized three-dimensional turbulence closure model which can be used 
in the calculation of practical engineering problems. This model overcomes a number 
of limitations of existing two-equation models. First, the current model can accu- 
rately predict the growth rates, velocity, and shear stress profiles for a variety of free 
shear flows using only one set of closure coefficients and boundary conditions. Fur- 
thermore, the fc-£ model allows for the solutions of a wall bounded flow and its wake 
without a modification of the model constants or boundary conditions. Moreover, 
the current model compares favorably to both the Johnson-King and k-ui models for 
the prediction of two-dimensional airfoils. For these airfoils, the current model does a 
good job predicting skin friction and shock location, as well as the extent and location 
of the separated regions. Also, the current form of the A;-£ model is free of damping 
functions and geometrical factors. With this development, one of the major obstacles 
in applying the model to three-dimensional flows has been removed. 

Because of the nature of the development process the current k-Q model still 
has one primary limitation. This is the models tendency to overpredict the eddy 
viscosity around a strong shock wave. The reason for this limitation is due to the 
fact that the pressure gradient term was calibrated for low speed flows (subsonic and 
transonic airfoils) with relatively weak shocks. When the model was used for flows 



51 


with strong shocks (Mach > 2.5) excessive values of fj. t were generated, and therefore, 
an overprediction of skin friction and heat transfer occur in the vicinity surrounding 
the shock. It is to be noted that, near the shock, over 50% of the turbulent energy 
production is a result of the work of the normal stresses. Since such stresses are 
not well predicted by two-equation turbulence models it is possible that the observed 
behavior of n t near a shock is a result of inaccuracies in the normal stresses. This 
conjecture can not be ascertained without a solution based on a stress model. 
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A Reynolds Stress Derivation 


In order to determine governing equations for the Reynolds stresses, moments are 
taken of the Navier-Stokes equations. This is accomplished by multiplying through 
by the fluctuating property and taking the time average of the product. 

The Navier-Stokes equations are re-written as 
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where jV is the Navier-Stokes operator. Take the following time average in order to 
derive an equation for the Reynolds stress tensor, 
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Consider each term separately, starting with the unsteady term 
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Next consider the convective term 
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and for the pressure gradient term 
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= 2/i 


<9x 


d / , , , , \ , du'i , du. 

9" ^mj o ^rr 


'dXr 


(A. 21) 


Combine Equations A. 9 - A. 21 to obtain the Reynolds Stress equation 


dt (Tl] ) + Urn dx 


dTjj _ dipUjUjU^) 

dx m 


dU x , , , dp TT ,dp 

Am r jm + +iq + U 3 

L/JLryfi L/ Jb j LS JL I 


+ 2/i 


d / , , / / \ / 


— s 


duj 

mi dIZ 


(A. 22) 


Rearrange the above equation into a more recognizable form. 


dr. 


dt 


V , rr 

i ^ m n 

uX m 


dr u ou 3 du, n ^ ac tjm 

+ Cij - fly + 


= An 


dXr 


- A 


jm 




dx„ 


n. = „■ + ?h. 

1] ySxj dxi 


( , du l , du. 


Cijrn — ""b P ~b P 'Uj&im 2/i ^U x S m j + UjS m ^J (A. 23) 


The above expression provides six new equations (one for each independent com- 
ponent of the Reynolds Stress tensor) which must be solved in order to completely 
determine the Reynolds stress tensor. 
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B Turbulent Kinetic Energy 
Derivation 

The equation describing the turbulent kinetic energy (per unit mass) is obtained 
by contracting the Reynolds Stress Eq.( A. 23). 

dC ilTn 


dr l: 

■ + U m 

dr n 

n dUi 


dt 

l dx m 

= 2 T lm 

dx m 

P^ii 

dk 

+ pUrr, 

dk 

dUi 

til 

dt 

1 3x 
L/u-m 

— ^im o 

dx m 

P J~ 



Dk 

dUi 

ptu 



P ~Dt 

7~im 0 
dXm 

2 


dxr, 


d 


dx m 

d 

dx m 


pU-U-Um 


pu'iUiU'rn 


+ P U m- 2 tiUiSn 


dk 


+ pu m - p 


0X T 


(B.l) 

(B.2) 

(B.3) 


By contracting the Reynolds stress equation we have reduced the tensor equation 
to a scalar one for r^. Since we are no longer dealing with a tensor equation, the 
scalar will be replaced with e. For an incompressible flow the pressure dilatation 
term, fl^, is zero and, therefore, if the Boussinesq approximation is used 


Tjj — 2 ptSij pkSij 


(B.4) 


the k equation can be written in its standard form for an incompressible fluid. 


dk TT dk 

P-ZI + PU, 


= TV, 


o , 1 t ' v m 1 nn 

dt dx m dx m 


dUi d 

- pt + 


dXr, 


dk 


PfaT ~ 2 P U 'i U 'i U 'm - P' u 'r 


(B.5) 


The above equation can be re-written in terms of the enstrophy as 


Dk Tij dUi d u t dk 

Dt p dxj dxj a* dxj J 


Xj 


(B.6) 
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From the enstrophv modeling we have 


= 


Vt 

2 G t 


d£l m 


dxi dx r 


+ 


L pmi 


d (u' p ufy 


dxi 


dk 

C)Xn 


(B.7) 


u'pu't = —2i/ t S p i + -k5 p i 


Therefore, after neglecting the higher order derivatives we obtain 


d_ 

dxj 



€pmi 9 

dk 

(«j 

6 dxj 

dx p 


(B.8) 


(B.9) 


Substitute the above to obtain the modeled k equation (in an incompressible 
medium) 


dk dk 

P-57 + pUm ~ — = Ti 


dt ** m dx im dx 

L/ 1 KJ Jb yji 


dUi d 

-K + 


dXr, 



P t_\ 

f 1 

[V3 


dx m 


(B.10) 


In order to relate the above equation to other two-equation models we can relate 
the dissipation e to the enstrophy, £. 


e = 2 us^ 

(B.ll) 

a 2 (u'«;) 
2 dxidxj 

(B. 12) 

= K+ 92 (“ iu 9 

S dxidxj 

(B.13) 


or for a large Reynolds number the above can be re-written as 1 


e = K 


(B.14) 
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C Vorticity Derivation 


Begin with the incompressible Navier-Stokes 


P 


dui 

dt 


+ Uj 


du j 
dxj 


dp dU j 

dxi dxj 


(C.l) 


where 


tij — 2/J.Sij 


(C.2) 


Re-write the above as 


dui ( dui 

P dt + U] dxi 


diii \ du, 

' 4 - 7 / - 

+ Uj dx t 


dp d 
dxi + P dxj 


' dui duj ^ duj 

dii 


( dui duj ' 
\dxj dxi t 


(C.3) 


or 


dui duj dp d 

+ 2r <n + + n—Clr,,) 


(C.4) 


Where is the rotation tensor and can be re-written as 1 


Tij = — ~€i jk^k 


(C.5) 


Now, substitute Equation C.5 into Equation C.4 to obtain 


du { duj 

^j^ijk^k “t* 


dp d 

+ tij k^h) 


dxi dx{ dxj 


(C.6) 


which can be re-written as 


dui d 

dt dii 


P UjUj 


~P + 2 


+ CijkUjUk - ue ijk 


duJk 

dxj 


(C.7) 
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From the definition of vorticitv 


U)i = 


^ ijk 


duk 

dxj 


(C.8) 


we see that in order to derive an equation for the vorticity you simply apply the curl 
operator to the above form of the Navier-Stokes (Eq. C.7). 


^ ilm 


d ( du. 


d 


dx m \ dt 


^ ilm 


dX T 


d 


_d_ 

dxi \p 

d u k 


P u j u j 


T ^ ilm 


d 


dx n 


( f-xjk UjUJk ) 


dx„ 


' ^ijk 


dxj 


(C.9) 


Note that the first term on the right hand side is zero from the multiplication of 
d2(£.+lPLL) 

a symmetric tensor ( — dXrndx \ - ) and a skew-symmetric tensor (e irm ). Therefore, the 
above reduces to 


duii 


dt dx 


d(ujUk) d 


dXri 


( duJk 

\ dxj 


(C.10) 


Apply the following identity 




mj 


(C.ll) 


to obtain 


_ djuiUrn) _ d^nUJi) _ _d_ ( dum \ _d_ / du t \ 

dt dx m dx m U dx m V dxi ) V dx m \ dx m ) 


The above equation can be simplified by noting that the divergence of the curl (V x V) 


is zero 


dut n 
dx„ 


= 0 


(C.13) 
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and by switching the dummy indices to obtain the instantaneous incompressible vor- 
ticity equation. 


dui z diJi dui d 2 u>i 

dt +Uj dxj Ll>3 dxj + U dxjdxj 


(C.14) 


or 


Dui 


dui d 2 u>i 

3 dxj dxjdx-j 


Dt 


(C.15) 
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D Enstrophy Derivation 

Start with the instantaneous vorticity equation for an incompressible flow (see 
Appendix C). 


Du)i dui d 2 u l 


ujj— h y- 


(D.l) 


Dt 3 dxj dxjdxj 
The first step is to derive the fluctuating vorticity equation. This is done by de- 
composing the instantaneous variables into a mean and a fluctuating component and 
then applying a time average. This will provide an equation for the mean vorticity. 
If you then subtract the mean equation from the instantaneous, you will obtain the 


fluctuating vorticity equation. 

Ut(xj, t) — Ui(x{) 

Substitute Equations D.2- D.3 into Equation D.l and apply a time average. 


(D.2) 

(D.3) 


dQ t , TT dn t , ,du[ „ dUi , ,du\ , d 2 Qi 

+Ui d^ +Ui d7 i =ni d^ +Uli a7 j + 


dt 


DQi 

Dt 


Qj 


dUi d 2 n { . du'i 
dxj 


+ y 




dui' 


- u\ 


dx) ' ] d Xj 3 dxj 


„ 3 u, . a 2 a . a r — , 

When the above equation is combined with the identity 

l] dx j ~ 3 13 


(D.4) 

(D.5) 

(D-6) 


(D.7) 
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the Mean Vorticity equation for an incompressible medium is obtained. 


DVLi d r -i—f -r-r 

~Di = " UiU > 


q , a 2 fi ( 

+ QjS t j + u q^2 


(D.8) 


Multiply Eq. D.8 by fli to obtain an equation for (mean equation) 

T-J Q ^2 Q 

— (DA) = -2 Sli— [^ - ^] + 2Q t Q J S l] + 2i/Q t 1 


dxj 


Dt l dx 3 

Multiply Eq. D.l by to obtain an equation for a;, a;, (instantaneous equation). 


(D.9) 


D , s 0 dui , 0 d 2 u)i 
— (a Wi) - 2^ l u j — + 2uu)i-^ 


(D.10) 


Dt 

Now decompose the above equation using Equation D.2- D.3 and apply a time average 
to obtain 


| + 1 (=R) = an^f + 2 ^g + + a^g + 2 ^g 


' dxj 


' dxj 1 J 


+2r/n^ + 2^'^ - ^ [aQi + 2D^'+ W 'a;']j 


o 

— Uj — — [Qi^i + + cj'cjJ] 

OXj 


(D.ll) 


Now subtract the mean and instantaneous equations to obtain a relation for cj'cj,' 


£(uM) = - 2 “W;g - £- (“»0 + + 2w;u,;s y + 2n,u ,'s' 


y 


a 2 


do;' 3a;' 


+1/ (ai'a;') — 2u- 

dx 2 m [ l l) dxj dxj 


(D.12) 


The above is the incompressible enstrophy equation where the enstrophy 


C = 


(D.13) 


is the sum of the squares of the fluctuating vorticity components. 
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E Modeling of the ( Equation 


We begin with the incompressible Enstrophy equation (see Appendix D) 


§ + ~ ~ 2u '^ ax, 


— da d ( u 'M^ 


dxj 


2u j'iUjS'ij + 2uj' i bj'jSij + 2 u}' i s' i jQ.j 


+v 


a 2 c 


2u 


du)[ dul 


dx j dxj dxj dx-j 


(E.1) 


From Eq. E.l we have six terms which require modeling. 


E.l Modeling of u[u' 


The first term is a second order tensor and is therefore re-written as the sum of a 
symmetric and skew-symmetric tensor. 

u'iU'j = \ + w'-uj) + ^ (uj'iu'j - (E.2) 

= Aij + B l} (E.3) 

E.1.1 Skew-Symmetric Portion of 


The first step is to examine the Skew-Symmetric portion of the above term. This 
is done by applying the following identity to Eq. E.3. 


d(u' m u') dk 

dxi dx m 

(E.4) 

~ tmij-A-ij + 

(E.5) 

“ 0 + 

(E.6) 
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Solve the above for B \j. 


tmii d _ {<X) dk 
2 dx t dx m 


(E.7) 


E.1.2 Symmetric Portion of u;' u'j 

The next step is to consider the symmetric portion of a;' u' . 

A ij = \ + u j u 'i) ( E - 8 ) 

(E.9) 
(E.10) 

= Uj ^ imt ^Im. ( E ‘^) 


ox m J 

, Ujml ( du\ dll' m \ ejrnl ( du\ g<,V 

Uj 2 \cto m dxi ) 2 \dx m dxi J 


— u j£iml r l m 


(E.12) 


u'jr[ m is a third order tensor. Therefore, in order to maintain the correct tensorial 
notation this is modeled using the gradient diffusion assumption. 


- v t dflim 
U ° Tlm ~ a T dxj 


Combine Eqs. E.7, and E.13 to obtain 

-J—J _ v^_ r (KU ernii d ( u 'm u 'i) _ dk 

UjU}l 2 a T dxj dxi 2 dxi dx m 


(E.14) 
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E.2 Modeling of 


du'-Ul'-LO 1 - 

111 


z'w'u;' is modeled using a gradient diffusion model. 




dxj 


-vt dC 

<7( dXj 

d_ Vt d( 

dxj dxj 


(E.15) 


(E.16) 


E.3 Modeling of wJa/sL — v 


T~r dwj du[ 
u dxjdxj 


u liUijS'ij - v is modeled using the suggestions given in References [47, 38] 


UiU'jS'ij - v 


dxj dxj t * 

iik) 


(E.17) 


where r is the turbulent time scale and R t is the turbulent Reynolds number ^ 7 . 
The above term was developed for high turbulent Reynolds number and homogeneous 
turbulence. Therefore, an additional term was added (for tu'w's'^) to account for low 
Reynolds number and non-homogeneous effects. 


wjwj-sy - v 


dui'i du[ /3 7 C 


t i _ 1 ^ n c ^ 

to, “ Q* ^75 


(E. IS) 


Based on the research to date, the term that gave the most satisfactory results for 
the airfoils under consideration was a result of vortex stretching. Because of this, the 
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above term, which was modeled assuming zero pressure gradient, in now modeled as 


lA-'V Ml , g >(»a) (“) R ‘ fE 

■ ‘ « ax, ax, - w n ' a ’ s ” SK + (• a!M | y (E ' 19) 

1 + g p v si ) 


E.4 Modeling of uj [ uj \ 


ujiijjj is modeled by noting the appropriate tensorial notation. Also, when con- 
tracted the right hand side must be equal to Q (from definition of £) 


uJiUJj = a 3 C b tj + -j- 


(E.20) 


_ -Ku'j 2 

bii ~ k + z 6i ‘ 


(E.21) 


where bij is the anisotropy tensor. 


E.5 Modeling of u-s-. 


ui'ts'ij is modeled by noting the appropriate tensorial notation. 


TTTT-s'.,. du 'm _ \ e Hm ( du' m du\ \ . € tlm ( C>< 

13 13 llm dxi 13 2 l dxt dx m ) 2 \ dx, 


T (E.22) 


— s ij ( e t ImS'mi + £Um r 'ml) 


(E.23) 


— s i 


(E.24) 


The above tensor, s' , r' ml , is symmetric in ’ij’ and skew-symmetric in ’ml’, therefore, 
is modeled as 


-r-^ = _ZP foCOnt ft ( 9k dC, Mj%_\ u t 
ij ml pk Q Q 2 V dx m dx t dx t dx m ) + ^v QQml (E>25) 
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where. 


Q = yj 


(E.26) 


. .i q i — ZlL 
^ " pk 


AC fit Ct/mA ( dk dC, 


Q. 


Q 2 V dxi dx r 


g-fj-fiSoD, 

OX m OXi I V 


(E.27) 


Combine all of the above terms to obtain the incompressible Enstrophy equation. 

d («{»«i) dk 
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(E.28) 


R k = 


Vy/C 


(E.29) 


This equation along with the Turbulent Kinetic Energy equation (see Appendix B) 

(E.30) 
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dx. 
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+ BL) 
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dx m 
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Ok) 

dx m 


comprises the current two equation model with the eddy viscosity calculated as 

(E.31) 
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F Free Shear Layers 


In order to solve the similarity form of the governing equations the momentum, 
turbulent kinetic energy, and enstrophy equations must be solved simultaneously. 
The momentum and turbulent kinetic energy equations are derived in Ref. [2] and 
are reproduced here for convenience. 


v dU _ ]_d_ 
dr] qJ dr] 


rfN— 

dr] 
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dK 1 d 


T? 


, N dK 


Ok dr] 


= S u U 


S k K + N 


dr] dr] 

where N is the non-dimensional eddy viscosity given by 

K 2 (i 7 ) 


'dLC 
dr] J 


- E 


(F.l) 

(F.2) 


N = a 


E{v) 


(F.3) 


and K(r /) and E(r]) are the non-dimensional turbulent kinetic energy and enstrophy, 
respectively. S u , S k , and V are provided in Table 4.1 Ref. [2]. The solution procedure 
used is that of Ref. [2] with the exception of replacing the e equation with the ( 
equation. Therefore, since the momentum and turbulent kinetic energy equations 
remain unchanged, the remainder of this appendix deals with the derivation of the 
similarity form of the equation for the wakes, jets, and mixing-layers. 

We begin with the enstrophy equation for a planar two-dimensional incompressible 
flow and by using the standard boundary layer approximation or thin layer approach. 


d( d( d 

U dx + V dy dy 
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d 2 u du du, 


dy 2 J dy 2 dy dy 
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F.l Plane Wake 
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u(x,y) = U , oo - 

k ( x ’V) = £#(»?) 


_ ai£!M ( p \ 

t_ Vpt/ooJ 


C(x,y) = ^^ 


Bin) 


N = C 

JV ^ £7(t7) 


n 


= y\[ 


pML 

Dx 


The first step is to non-dimensionalize the equations (using the above) and re-write 
the partial differential equation as a total differential equation. 


d 

iNdE] 

„ & E 2 „ r „ , 
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+ 3&JV 


dU 


di 7 


(F.5) 


The next step is to solve the ODE using the method described in Reference [2]. 
This includes applying the Rubel-Melnik transformation 48 to reduce the stiffness of 
the equation set. This stiffness can be explained by noting that at the edge of the tur- 
bulent region the equations become stiff because both k and £ approach zero. Since we 
are not interested in solving the equations at the turbulent/non-turbulent interface, 
we will shift this boundary to infinity by applying the Rubel-Melnik transformation. 


d d 
rff ''‘dy 


(F.6) 
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or. in terms of the non-dimensional eddy viscosity 


d£ dr] 


From the above transformation, it is seen that as the eddy-viscosity approaches zero, 
the <f coordinate will approach infinity, therefore, we simply integrate the equa- 
tions from the centerline out to some prespecified f max which occurs prior to the 
turbulent/non-turbulent interface. 

After applying the Rubel-Melnik transformation the equation for the plane wake 
takes the form 


,dE d f 1 dE 


ft E 2 N a 3 CuK (dU\ 2 2 


406 dl£ 3 
+ 3iV ~dt 

A similar procedure is used for the Mixing-Layer and the jets. 


F.2 Mixing Layer 


u (x,y) — UiU(t]) v(x, y) = Ui [V + rjU\ 
k(x,y) = U 2 K(rj) ((x,y) = E(r ] ) 
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*1«S 
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After the Rubel-Melnik transformation the above equation becomes 
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F.3 Jets (Plane, Round, and Radial) 

The two-dimensional axisvmmetric form of the enstrophy equation is 
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where the non-dimensionalization parameters are 

u(x, y) = U 0 U (t}) v(x , y) = U a V (77) 


U(rj) = 


_ {jTIl 


V(T}) = TjU- F 2 j - 1 
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(F.10) 


(F 11) 


j = 0, m = 0 Planar Jets 


(F.12) 
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7 = 1, m = 0 Radial Jets 
j=l, m = 1 Round Jets 


(F.13) 

(F.14) 


Substitute in the above and re-write in Dr. Wilcox’s notation (V = —2 to 
obtain 
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After the Rubel-Melnik transformation the above equation becomes 
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G Log Layer Analysis 

The total shear stress is represented as the sum of the laminar and turbulent 
stresses. 

r = T t +T t (G.l) 

From experimental measurements it has been determined that the laminar shear stress 
accounts for approximately 5% of the total stresses in the log law region. Therefore, 
in this region we will assume that the laminar shearing stresses are negligible when 
compared to the turbulent shearing stresses. It is also noted that the total stress is 
nearly constant in the log law region. As a consequence of the above assumptions we 
will approximate the total shearing stress in the log law region as 

r = n (G.2) 

= r w (G.3) 

where t w represents the shear stress evaluated at the wall. Also from experiment, it 
has been determined that the turbulent kinetic energy, k, is approximately constant 
in the log law region. Therefore, we can represent the ratio of the turbulent shearing 
stress to the turbulent kinetic energy as 49 

~ = \Jc~n where , C M = 0.09 ( structural factor ) (G.4) 

In the log law region the velocity profile is given from the law of the wall, 

U + ~ — lny + + B 

K 


(G.5) 
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where k is the Karman constant (« .40) and B is a dimensionless constant (« 5.0). 
The velocity and normal distance in Eq. ( G.5) were non-dimensionalized using the 
friction velocity, u T , where 


e 

H 

II! 

^i?l 

(G.6) 

Therefore, 


u^ u - 
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(G.7) 

+ _ U rV 
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v 
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The Reynolds stress tensor in the log law region is approximated as 


T t —pUV 

(G.9) 
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(G.10) 

where the eddy viscosity, u t is given by 


Pt 

"t = = , 
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(G.ll) 

Combine the above relations to obtain 



(G.12) 
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k 
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Therefore 



Similarly, by solving for £ in the log law region we obtain 



vKy 


and the eddy viscosity becomes 


v t = u T Ky 


(G.16) 


(G.17) 


(G.18) 


Next we will examine the vorticity for a two-dimensional cartesian coordinate system. 

a, = (G.19) 

For two-dimensional flows the above reduces to 


By also noting that 


Q x = Q 2 = 0 
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(G.24) 


Similarly, the mean strain rate is approximated as 
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1 dU 
. 

2 dy 
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2/cy 

and the component of the rotation tensor becomes 


(G.26) 

(G.27) 
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2ny 

Substitute the above relations into the Enstrophy equation in order to obtain a rela- 
tion between the constants. 


(G.28) 

(G.29) 


-UjU, 


7 — rdQi —u 3 1 


dxj 2 k<j t y 3 


d /, i ,\ ut 1 


- V 


v<*c, y 

du)\du 

dxj dxj Ai/K?y 2 


— 7 — T c a 3 U \\[Cp- 




z/ 


a 2 C = 

dxjdxj k y 3 


Combine the above to obtain 

]_ f u t /?5 a$u* s J(T tl u *p 4 2/3 6 il A 

y 2 \ 2ua q Auk 2 Auk. 2 3uk 2 3i/k 2 I 


(G.30) 

(G.31) 

(G.32) 

(G.33) 

(G.34) 

(G.35) 


+hot = 0 


(G.36) 
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From a power series analysis we note that each of the above coefficients must be zero. 
The p- term is considered a higher order term and is therefore neglected. Therefore, 
the log-law analysis has provided us the following relation for the model constants. 

Qk 2 

^ ~ [zihsfC^ - 3*^ - 4(3 4 - 8ft] 


(G.37) 
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H Boundary Layer Analysis 

Begin with the Reynolds averaged equations of motion for an incompressible, two- 
dimensional, planar flow 


dxi 

dUj dU^ _ _dP_ dt x± 
^ dt + ^ 3 dxj dx t dx.j 


(H.l) 

(H-2) 


then by making the boundary layer assumptions the above equations reduce to the 
standard form of the boundary layer equations 


dP d \. .dU 


(H.3) 

(H.4) 


This equation set is parabolic in the x direction and, therefore, is solved using a 
marching procedure as described in Reference [15]. This reduction in the complexity 
of the equation set reduces the typical run times from several hours to approximately 
one minute. This reduction in computational effort allows a wide variety of modeling 
and closure coefficients to be tested. 


The model was added to the Harris and Blanchard 15 boundary layer code 
in place of an existing k-ui turbulence model. The only difference in the solution 
procedure used in the current research as compared to the original code 15 is seen 
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from the non-dimensionalization of which uses 


C = 


Nr = 


UrNrC 

L r 

UrLrPr^ 

Hr 


(H.5) 

(H.6) 


where R, *, N ^ represent a reference quantity, a non-dimensional quantity, and a 
reference Reynolds number, respectively. The reference Reynolds number is used 
to scale the Enstrophy in order to maintain numerical accuracy during the solution 
procedure. 
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I Near Wall Behavior of k and £ 

Consider the turbulence quantities as you approach a solid surface by expanding 
the fluctuating velocity components in a power series 

u '{y) — a o + a \y + a 2 y 2 + a 3j/ 3 + ••• 

v'(y) = b 0 + b x y + b 2 y 2 + b 3 y 3 + ... (1.1) 

w'{y) = c 0 + c x y + c 2 y 2 + c 3 y 3 + ... 


Now impose the conservation of mass for an incompressible fluid 

du L = d<= o 

dxi dx t 

and the ’no slip’ boundary condition 


(1.2) 


u'(0) = v'(0) = w'( 0) - 0 (1.3) 

to obtain 


u'(y) = a x y + a 2 y 2 + a 3 y 3 + ... 

v'(y ) = b 2 y 2 4- b 3 y 3 + ... (1.4) 

™'{y) = C\y + c 2 y 2 + c 3 y 3 + ... 

Therefore, from the definition of turbulent kinetic energy we obtain 

k = = \ (a\y 2 + bly A + c 2 y 2 ) 

- k 0 y 2 as y — ► 0 


(1.5) 

( 1 . 6 ) 
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where k 0 is some constant. 

A similar analysis is performed for the enstrophy at a solid surface. 


c = <44 

_ ” du 'k ~ 

a *ilm c\ 

OXj OX i 

( du' \ 2 / dw' \ 2 

V By ) + V dy ) 

= + c]) + ... 

= Co asy -> 0 


(1.7) 

(1.8) 

(1.9) 


( 1 . 10 ) 


where Co is some constant. 
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J Homogeneous shear flows 


Sij = 


A homogeneous shear flow 

r 

0 5 0 
5 0 0 
0 0 0 

is given by enforcing a linear velocity profile with no dependence on the x or z coor- 
dinate directions. In other words, the following must be imposed 


dU „ 

— — = Constant 
dy 

(j.i) 

dU _ dU _ o 
dx dz 

(J.2) 

o 

II 

£ 

cb 

(J.3) 


where U is the mean velocity. 

The above constraints reduce the k and C, equations to 


dk _ ( du 

~di~ Ut \d^ 

_ c*3 C^k 
dt v 



&C 2 \/2 

(^+i) 




du 


du 

dy 

3 u 2 C, 

dy 


(J.4) 

(J.5) 


Now apply the following non-dimensionalization (* superscript represent a non- 


dimensional quantity) 


k — k 0 k 


(J.6) 
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CSk 0 
2 v 


f 



to obtain 


(J.7) 

(J.8) 

(J.9) 


dk * 
dt* 

dc 

dt‘ 


*)2 


C„(h-) 


C* 

= q 3 c„/f - 


-c 


ACC *) 2 


, , 4ftC„((c-) 2 

3 3 <• 


(J.10) 

(J.ll) 


The above ODE’s are solved using a fourth order runge-kutta integration scheme 
to determine the decay rate of turbulence for a homogeneous shear flow. An initial 
condition , corresponding to setting production equal to dissipation, gives a homoge- 
neous shear parameter of at t = 0, to begin the integration. 

It should be noted that many modelers use this type of flow to calibrate closure 
coefficients, however, since we were primarily interested in wall bounded shear flows 
and free shear layers we chose not to use the decay of homogeneous turbulence as a 
method of determining constants. Instead, we used homogeneous shear flow as a test 
case for the current two-equation model. 

The data from our test case includes an examination of the shear parameter > 
anisotropy tensor ( bij ), and the predicted decay rates of k and C 35 P er suggestion of 
Reference [37]. 
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K Free Stream Boundary 
Condition 


In this case, the equation governing k and £ can be written as 


dk 

dt 

dC 

dt 


31 = -K 


= -A 


(r?+s) 


(K.l) 

(K.2) 


where 


_ 1 ^ 
ft = “ 27 = 73- 

i/ 2 C C> i/ 


(K-3) 


Dividing Eq. (K.l) by Eq. (K.2) gives 


dA: A: + <5i/ 


d( &C fcC* 

Equation (K.4) is a first order linear equation. Its solution can be written as 


(K.4) 


d 


* = + C 


(K.5) 


(A " 2) 

where C is a constant to be determined from initial conditions. It can be seen from 


Eq. (K.5) that 


k 


25 


£r (*"0 


+ -C 


(K.6) 


(/?5 — 2 ) 1 / 

For p 5 = 2.37, the above equations show that as £ -* 0, and v t increase without 
bound. Therefore 


C = 0 


(K.7) 
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and initial values of £ must be chosen such that 

7 =^^ c -((^ 2)) 2 < K - 8 > 

Thus, for decaying homogeneous turbulence k and £ approach zero in such a way that 
^ remains constant. This, in turn, implies that during the decay process 

£ oc k 2 (K.9) 

Substituting Eq. (K.9) into (K.l) yields 

k oc t~ l , £ oc t ~ 2 (K.10) 
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L Final k-( Equation Set 


The final form for the modeled k - £ equations can be written as 
Dk dU, , d 


^ Dt im dx m ^ + dx m 

DC, p t dQi ( dQi <9fL 

p— = — I - — + 



Pt} 

dk 

u + 

VkS 

dx m 


r n k 
— Ci p — 


Dt a r dxj \ dxj 

Pfh 


dxi 


cqpvm | d 

Tp dxj 


( ,Pt\dC 
P + — 


a £ J dxj 


Rk + 6 

+ max 


, f ru , 2 r c MmMj 

C 2 + [azpCbij + -bijpc j Sij — 

,o.ol 





2 ^ 'Minn-n + 

kv Q? 


+ max 


djpUjP) 

dx t 


, 0.0 


kD 


J v P (Jp (1 + dp) 


(L-l) 


(L.2) 


where 


T{j — 2 pi 


Sij - fs kk 


2 r i r ^ 

-~5 tJ pk , v> t = C^— 


S 2 = S ij S ij , Q 2 = a^« , Rk = 


Vy/C, 


, Rt — R k 


5p = —a 

P 


kR t ( dp 


\| C 


, Cx = 0.60 , C Cl = 2.10 


_1_ _ 1 
a ~ p\ 



k 

Uf 
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9 Tables 


Table 9.1: Free Shear Flow Spreading Rates 


Flow 

Enstrophy 

Epsilon 

Measured 

Far Wake 

0.3130 

0.256 

0.365 

Mixing Layer 

0.1054 

0.098 

0.115 

Plane Jet 

0.1143 

0.109 

0.100 - 0.110 

Round Jet 

0.0906 

0.120 

0.086 - 0.095 

Radial Jet 

0.0965 

0.094 

0.096 - 0.110 


Table 9.2: Turbulence Model Closure Coefficients - Original 


Constants 

k- C 

<*3 

0.35 

P* 

0.42 

A 

2.37 

A 

0.10 

$7 

1.50 

As 

1.15 


0.6 

<7p 

5.6 


1.80 

x 

1.46 

c ul 

32.0 

C U 2 

32.0 

5 

0.10 
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Table 9.3: Turbulence Model Closure Coefficients - Final 


Constants 

k- C 

«3 

0.35 

A 

0.42 

05 

2.37 

06 


07 


08 

1.15 

O r 

0.07 

o v 

0.065 

o p 

91.9 

IHEH 


X , 


5 



Table 9.4: Cases Considered 


Flat Plate 

Mach # 

Re 


Grid 


0.20 

9.0x10* 


NAx201 (BL code) 

Airfoil 

Mach # 

Re 

AOA 

Grid 

NACA 0012 

0.50 

2.91xl0 b 

2.06° 

287x81 

| NACA 0012 

0.799 

g.ozio 6 

0.00° 

287x81 

RAE 2822 

0.73 

6.5a: 10 6 

2.80° 

321x91 

j RAE 2822 

0.73 

6.5xl0 6 

2.80° 

231x51 

RAE 2822 

0.73 

6.5xl0 6 

2.80° 

143x41 j 

NACA 4412 

0.20 

1.5xl0 e 1 

13.87° 

237x91 

NACA 0012 

0.55 

9.0x10* 

8.34° 

287x71 

NACA 0012 

0.799 

9.0x10* 

2.26° 

287x71 

RAE 2822 (case 9) 

0.73 

6.5x10* 

2.80° 

321x91 

RAE 2822 (caselO) 

0.75 

6.2x10* 

2.72° 

231x61 

Cylinder-Flare 

Mach # 

Re 


Grid 


2.89 

15.0x10* 


87x99x99 








10 Figures 



Figure 10.1: Wake flow schematic 






Figure 10.4: Self-similar turbulent kinetic energy (Planar Wake) 
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Figure 10.5: Schematic diagram of planar jet 






















Figure 10.19: Near wall turbulent kinetic energy for flat plate 



Figure 10.20: Near wall Enstrophy for flat plate 


104 



Figure 10.21: Near wall shearing stress for flat plate 



Figure 10.22: Flat plate log-law constant 










Figure 10.26: Near wake growth rate 
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y/b 

Figure 10.30: Wake shear stress profiles 



Figure 10.31: Time decay of Homogeneous shear parameter 



Figure 10.33: " 
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Fieure 10.34: Time decav of c + 



Figure 10.35: Boundary condition comparison of pressure coefficient for NACA 
0012 (Moo = -502, Roo = 2.91e6, a = 2.06) 


K (Enstrophy, bcl > 
K (Enstrophy, bc2) 
k -to (Wilcox} 


Figure 10.36: Boundary condition comparison of skin friction coefficient for NACA 
0012 (Moo = .502,^ = 2.91e6,a = 2.06) 


k< (Enstropfty, Del ) 
kn; (Enstrophy, bc2) 
k-« (Wilcox) 


Figure 10.37: Boundary condition comparison of pressure coefficient for NACA 
0012 (Moo = 0.8, Roo = 9.0e6, a = 0.0) 






114 


321x91 

231x51 

143x41 

□ Experimental 



i 

0.8 


RAE 2822 grid study 




Figure 10.41: Airfoil geometry for NACA 4412 
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Figure 10.42: NACA 4412 Pressure coefficient 



Figure 10.43: NACA 4412 Skin friction coefficient 
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Figure 10.44: NACA 4412 Velocity profiles 





x/c 


Figure 10.46: NACA 0012 (Case 1) Pressure coefficient 



Figure 10.47: NACA 0012 (Case 1) Skin friction coefficient 
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Figure 10.48: NACA 0012 (Case 2) Pressure coefficient 



Figure 10.49: NACA 0012 (Case 2) Skin friction coefficient 
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Figure 10.54: RAE2822 (Case 10) Skin friction coefficient 
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Figure 10.55: Geometry of cylinder /offset-flare juncture 
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Figure 10.58: Streamwise pressure comparison for 9 = 0 plane 



Figure 10.59: Streamwise pressure comparison for 6 = 90 plane 
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Figure 10.62: Skin friction coefficient for 8 = 90 plane 



Figure 10.63: Skin friction coefficient for 8 = 180 plane 



